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The  purpose  of  this  lecture  is  to  assess  the  current  state  of  our  ability  to  simulate  the  consequences  of  a 
fire  in  a  large  building,  and  suggest  some  areas  where  improvement  is  needed.  Attention  is  focused  on  the 
coupling  of  fire  dynamics  simulations  and  heat  transfer  analyses  to  each  other  and  to  structural  analyses  of 
the  damaged  building.  The  role  that  uncertainty  in  “input  parameters”  resulting  from  coupling  a  sequence 
of  complex  simulations  is  considered.  The  methodology  used  in  the  NIST  investigation  into  the  collapse 
of  the  World  Trade  Center  Towers  will  be  described  from  this  perspective.  The  intent  is  not  to  summarize 
the  results  of  the  investigation,  but  rather  to  provide  a  specific  context  that  illustrates  the  strengths  and 
weaknesses  of  the  methodologies  employed.  Research  needs  are  emphasized  by  examination  of  some 
basic  problems  in  fire-structure  interactions. 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

There  has  been  a  resurgence  of  interest  in  the  response  of  build¬ 
ing  structures  to  fires  over  the  past  several  years.  This  interest  was 
greatly  enhanced  by  the  attack  on,  and  subsequent  collapse  of,  the 
World  Trade  Center  ( WTC)  towers.  This  has  resulted  in  at  least  three 
detailed  physics-based  analyses  of  the  collapse  of  the  towers  in  the 
literature.  Usmani  et  al.  (2003)  and  Hori  (2004)  concentrated  on 
the  degradation  of  the  towers  given  assumed  temperature  distribu¬ 
tions  in  the  load  bearing  structures.  Abboud  et  al.  (2003)  performed 
a  more  comprehensive  analysis,  dealing  with  both  the  initial  impact 
of  the  aircraft  as  well  as  the  final  collapse,  using  as  much  photo¬ 
graphic  and  video  evidence  as  was  available  at  the  time  (Beyler 
et  al.,  2003 ;  Thater  et  al.,  2003 ).  More  studies  are  surely  on  the  way. 
NIST  has  just  completed  its  own  investigation,  the  details  of  which 
are  available  at  the  website  http://www.wtc.nist.gov. 

Rather  than  debate  the  causes  of  the  collapse  of  the  towers,  it 
seems  more  useful  at  present  to  discuss  the  methodology  used  to 
study  the  collapse.  Many  if  not  most  of  the  issues  that  had  to  be 
addressed  in  the  WTC  simulations  would  arise  in  the  analysis  of 
any  major  fire  affecting  the  structural  stability  of  a  large  building. 
The  complexity  of  such  buildings,  together  with  the  need  to  couple 
phenomena  usually  considered  in  isolation,  forces  compromises  in 
methodology  not  required  for  simpler  fire  scenarios.  This  statement 
certainly  applies  to  the  analyses  cited  above,  as  well  as  the  NIST 
investigation.  Hopefully,  the  discussion  will  shed  some  new  light 
on  our  ability  to  simulate  such  events. 
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First,  the  sequence  of  simulations  of  the  consequences  of  the 
attacks  on  the  towers  and  the  software  employed  are  discussed 
briefly.  Very  few  details  will  be  given,  as  the  purpose  is  primar¬ 
ily  to  orient  the  reader.  Since  the  software  employed  is  widely 
available  either  as  commercial  products  or  NIST  freeware,  it  is 
fair  to  assume  it  is  representative  of  the  current  state  of  the  art 
(indeed,  it  was  chosen  primarily  for  this  reason).  This  is  followed 
by  a  discussion  of  the  approximations  needed  to  couple  the  vari¬ 
ous  elements  of  the  simulation  together.  The  role  of  fire-induced 
uncertainty  and  the  extent  to  which  the  descriptions  of  the  different 
physical  processes  must  be  coupled  will  be  emphasized.  A  sam¬ 
pling  of  results  (drawn  primarily  from  the  cited  references)  will  be 
used  to  illustrate  these  points.  The  paper  will  conclude  with  some 
recommendations  for  future  research.  This  work  is  an  extended 
update  of  the  Emmons  Plenary  Lecture  presented  by  the  author  at 
the  Eighth  International  Symposium  on  Fire  Safety  Science  (Baum, 
2005). 


2.  The  NIST  approach— a  synopsis 

The  NIST-led  analysis  of  the  events  initiated  by  the  attack  on 
the  towers  is  composed  of  four  parts;  the  initial  impact,  the  fire 
dynamics,  thermal  analysis  of  the  load-bearing  structure,  and  struc¬ 
tural  deterioration.  The  information  flow  in  principle  proceeds  as 
follows:  The  initial  impact  simulation  is  needed  as  input  into  the 
following  three  parts.  It  partially  defines  the  geometry  used  for  the 
fire  dynamics  simulation,  provides  guidance  about  insulation  dam¬ 
age  that  affects  the  heat  transfer  calculation,  and  determines  the 
state  of  the  structure  that  survived  the  impact.  The  fire  dynam¬ 
ics  provides  the  thermal  environment  in  the  gas  phase  needed  to 
determine  the  heat  transfer  to  the  exposed  building  surfaces.  The 
thermal  analysis  determines  the  temperature  distribution  in  the 
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load  bearing  structure.  Finally,  the  structural  analysis  integrates  all 
the  previous  information  to  make  predictions  about  the  loads  and 
deflections  up  to  the  point  of  collapse.  The  structural  dynamics  dur¬ 
ing  the  final  collapse  of  the  towers  is  not  part  of  the  investigation. 

2.2.  Impact  analysis 

Simulation  of  the  initial  impact  of  the  aircraft  on  the  towers, 
performed  under  contract,  covers  approximately  the  first  0.6  s  of 
elapsed  time.  The  simulations  were  performed  using  the  com¬ 
mercial  code  LS-DYNA,  “a  general  purpose  finite  element  code  for 
analyzing  the  large  deformation  dynamic  response  of  structures 
including  structures  coupled  to  fluids”  (Hallquist,  1998).  It  should 
be  noted  that  although  the  code  has  some  fluid  mechanics  simula¬ 
tion  capability,  it  cannot  simulate  combustion  phenomena.  Thus, 
no  attempt  was  made  to  simulate  the  initiation  of  the  fireball 
in  this  part  of  the  investigation.  LS-DYNA  and  the  finite  element 
code  FLEX  (Vaughan,  2002)  used  by  Abboud  et  al.  share  a  common 
ancestor  DYNA-3D,  although  both  are  undoubtedly  much  more 
capable  now.  The  codes  both  used  highly  detailed  representations 
of  the  aircraft  structure  and  those  floors  of  the  WTC  towers  in  the 
immediate  vicinity  of  the  impact  to  produce  an  estimate  of  the  dam¬ 
age  inflicted  on  the  load  bearing  structure  by  the  breakup  of  the 
aircraft. 

The  predicted  damage  to  the  exterior  of  the  towers  was  veri¬ 
fied  by  direct  comparison  with  the  highly  detailed  photographic 
and  video  evidence.  The  initial  configuration,  material  properties, 
and  relative  velocity  of  both  aircraft  and  building  are  known  quite 
accurately.  The  level  of  sophistication  in  the  computer  codes  then 
permits  an  impressively  detailed  prediction  of  the  breach  of  the 
outer  perimeter  of  the  towers.  However,  the  most  important  out¬ 
puts  of  the  analysis  are  the  prediction  of  the  post-impact  interior 
geometry  and  the  degree  of  damage  to  the  insulation.  For  the  pur¬ 
poses  of  the  present  paper,  it  completes  the  description  of  building 
needed  for  subsequent  analyses.  Thus,  the  remainder  of  this  work 
will  regard  the  post-impact  state  of  the  building  as  defined  by  these 
calculations  as  the  initial  conditions  for  the  fire  dynamics,  thermal, 
and  structural  analyses. 

2.2.  Fire  dynamics 

The  fires  induced  by  the  aircraft  impact  were  studied  using 
the  NIST  fire  dynamics  simulator  (FDS),  a  finite  difference  based 


computational  fluid  dynamics  code  solving  the  low  Mach  number 
form  of  the  Navier-Stokes  equations  appropriate  for  combustion, 
smoke  and  heat  transport  generated  by  fires  (McGrattan,  2004). 
FDS  was  modified  to  enable  parallel  processing  as  part  of  this  study 
(McGrattan  and  Bouldin,  2004).  The  code  uses  a  geometry  model 
of  those  floors  of  the  towers  affected  by  fires.  The  geometry  model 
differs  substantially  from  that  used  for  the  impact  analysis,  since 
it  is  based  on  architectural  rather  than  structural  features.  The 
spatial  resolution  employed  is  too  coarse  to  include  much  of  the 
load-bearing  structure,  retaining  only  the  floor  slabs.  The  core  and 
exterior  columns  are  present  but  cannot  be  represented  to  scale. 
On  the  other  hand,  it  includes  the  office  partitions  and  work  sta¬ 
tions.  The  model  also  includes  vents  opened  up  either  by  estimates 
of  damage  caused  by  the  initial  impact  or  windows  broken  dur¬ 
ing  the  fire.  Fig.  1  shows  the  geometry  used  by  FDS.  Every  feature 
shown  in  the  figure  is  included  in  the  simulations.  No  attempt 
was  made  to  predict  window  breaking.  Photographic  evidence  was 
used  to  obtain  times  and  locations  of  window  breakage  and  the 
resulting  information  was  specified  as  part  of  the  input  boundary 
conditions. 

The  fireballs  resulting  from  the  initial  impact  were  studied  sep¬ 
arately  to  estimate  the  fuel  consumed  outside  the  towers  and  thus 
unavailable  as  an  accelerant.  Two  independent  techniques  were 
used  to  obtain  estimates  of  both  the  fuel  consumption  and  the 
uncertainty  of  the  estimate.  FDS  simulations  were  performed  by 
Rehm  et  al.  (2002)  based  on  a  variety  of  assumptions  about  the 
internal  configuration  of  the  towers  following  impact.  Baum  and 
Rehm  (2005)  performed  an  analytical  study  assuming  the  shape  of 
the  fireballs  emerging  from  each  face  of  the  towers  to  be  hemi¬ 
spherical.  Comparison  of  these  calculations  with  video  images 
yielded  estimates  that  10-25%  of  the  aircraft  fuel  was  consumed. 

The  FDS  simulations  of  the  fires  for  the  remaining  life  of  the 
towers  were  tested  by  comparing  predicted  temperatures  and 
heat  release  rates  with  those  measured  in  a  series  of  labora¬ 
tory  experiments.  These  experiments  burned  paper  laden  work 
stations  similar  to  those  used  in  the  WTC  towers  in  a  compart¬ 
ment  with  limited  openings  resembling  those  in  the  towers.  The 
actual  office  contents,  their  condition  following  the  impact,  and 
the  completeness  of  burning  are  additional  sources  of  uncertainty. 
Many  simulations  of  the  actual  fires  in  the  towers  with  varied 
assumptions  about  the  input  parameters  were  performed  to  test 
the  sensitivity  of  the  model.  Calculated  flame  locations  near  the 
perimeter  were  compared  with  photographic  and  video  images. 


NET 


Fig.  1.  Layout  of  North  WTC  tower  model  as  used  in  FDS  simulations.  Windows,  interior  partitions,  and  workstation  geometry  are  all  included.  However,  the  columns  are 
not  to  scale  and  no  internal  floor  system  structural  components  are  present  in  the  simulations. 

Courtesy  McGrattan  (2004) 
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Fig.  2.  The  perimeter  column  structure  of  North  WTC  tower  shown  as  simulated  in 
the  thermal  analysis  using  ANSYS.  The  details  of  the  column  shapes  and  connection 
with  spandrels  are  shown  as  inserts. 

Courtesy  I<.  Prasad  (Baum,  2005). 


Fig.  3.  The  North  WTC  tower  floor  system  (with  concrete  slab  removed  for  clar¬ 
ity)  and  core  columns  shown  as  simulated  in  the  thermal  analysis  using  ANSYS. 
The  details  of  the  trusses  including  the  insulation  are  shown  as  inserts.  Additional 
illustrations  can  be  found  in  Prasad  and  Baum  (2005). 

Courtesy  K.  Prasad  (Baum,  2005). 


Details  of  these  and  other  aspects  of  the  WTC  fire  simulations  can 
be  found  in  McGrattan  and  Bouldin  (2004). 

2.3.  Thermal  analysis 

The  thermal  analysis  of  the  structure  couples  the  fire  dynam¬ 
ics  to  the  calculations  of  stress  and  deflection  that  ultimately 
determine  the  collapse  path  of  the  towers.  The  analysis  itself  is 
composed  of  two  parts;  the  calculation  of  the  local  heat  trans¬ 
fer  from  the  heated  gases  to  the  exposed  surfaces  of  insulation 
materials  and  structural  components,  and  the  heat  flow  through 
the  condensed  phase  materials.  A  decision  was  made  to  use  the 
code  ANSYS  (2003)  to  perform  the  structural  analysis.  Since  ANSYS 
also  has  a  thermal  analysis  capability,  this  led  the  development 
of  a  methodology  that  couples  ANSYS  to  FDS.  This  methodol¬ 
ogy,  the  fire-structure  interface  (FSI),  is  capable  of  generalization 
to  other  codes  and  other  structural  fire  problems.  It  makes  use 
of  a  mixture  of  analytical  and  computational  techniques  to  cope 
with  the  wide  disparity  of  length  and  time  scales  that  control  the 
physical  processes  simulated  in  FDS  and  whatever  structural  anal¬ 
ysis  code  is  chosen  for  use.  A  more  detailed  description  of  the 
coupling  procedures  that  comprise  FSI  is  presented  in  the  next  sec¬ 
tion.  These  techniques  are  also  summarized  in  Prasad  and  Baum 
(2005). 

The  finite  element  model  used  for  the  ANSYS  simulation  of  the 
heating  of  the  structure  is  the  most  detailed  of  all  those  used  for  the 
investigation.  It  is  composed  entirely  of  three-dimensional  “brick” 
elements,  and  represents  all  the  beams,  columns,  and  floor  systems 
on  all  floors  affected  by  the  fires  (including  the  floors  immediately 
above  and  below  floors  on  which  fires  occured).  It  incorporates 
detailed  models  of  the  protective  insulation,  including  estimates 
of  damage  obtained  from  calculated  debris  patterns  performed  as 
part  of  the  impact  analysis.  Random  irregularities  in  the  thickness 
of  spray-on  insulation  are  also  taken  into  account.  Fig.  2  shows  the 
perimeter  columns  and  connections  as  simulated  for  the  thermal 
analysis.  The  floor  system  and  core  columns  are  shown  in  Fig.  3. 
Analogous  models  were  made  for  the  South  tower.  The  accuracy 
of  the  calculations  was  tested  by  comparing  FSI  predictions  with 
a  series  of  large  scale  experiments  on  steel  trusses  and  columns. 
Studies  of  this  type  have  been  performed  by  Flasemi  and  his  col¬ 
laborators  (Hasemi,  2000).  The  results  indicate  that  FSI  can  reliably 
predict  the  temperature  in  structural  components  and  that  it  does 


not  add  any  more  uncertainty  to  the  overall  analysis  than  that 
already  present  in  the  FDS  simulations.  A  detailed  analysis  of  the 
experimental  uncertainty  and  model  sensitivity  associated  with 
both  the  fire  dynamics  and  thermal  analysis  of  the  laboratory  exper¬ 
iments  was  also  performed  (Hamins  et  al.,  2005). 

This  level  of  detail  is  necessary  precisely  because  the  calcula¬ 
tions  it  supports  provide  the  transition  from  large  scale  convective 
and  radiative  transport  in  the  gas  to  highly  localized  thermal  con¬ 
duction  through  solids.  It  is  also  needed  because  it  turns  out  that 
the  temperature  distributions  in  the  steel  are  more  sensitive  to  the  fire¬ 
proofing  thickness  and  the  degree  of  impact  damage  to  the  insulation 
than  any  other  quantity.  It  is  feasible  only  because  heat  transfer  by 
thermal  conduction  (even  non-linear  conduction  with  non-linear 
radiative  boundary  conditions)  is  much  simpler  computationally 
than  either  turbulent  combustion  or  non-linear  solid  mechanics. 

2.4.  Structural  analysis 

The  stresses  and  displacements  in  the  structure  were  calculated 
using  ANSYS  by  a  contractor  with  the  thermal  input  from  FSI  and  the 
initial  impact  damage  estimate  obtained  from  LS-DYNA.  A  series  of 
highly  detailed  component  models  of  relevant  parts  of  the  structure 
were  subjected  to  a  variety  of  loads.  These  were  used  to  develop  a 
less  detailed  set  of  floor-level  models  as  well  as  multi-floor  mod¬ 
els  of  portions  of  the  perimeter  column  system.  Video  images  were 
used  to  select  critical  portions  of  the  perimeter  columns  for  detailed 
analysis.  Finally,  a  global  model  based  on  simplifications  derived 
from  sensitivity  analyses  using  the  component  and  floor  level  mod¬ 
els  was  developed.  It  should  be  noted  that  local  models  and/or 
highly  simplified  floor  system  models  were  also  used  by  Usmani 
et  al.  (2003)  and  Hori  (2004)  in  their  analyses. 

An  interpolation  scheme  (part  of  FSI)  was  used  to  input  the  pre¬ 
dicted  temperatures  into  the  finite  element  models  employed  for 
the  structural  analysis.  The  thermal  information  was  transferred  to 
the  ANSYS  structural  model  at  intervals  of  1 0  min  of  simulated  time 
for  each  of  the  towers.  Neither  Usmani  nor  Hori  attempted  this  step. 
This  process  introduces  additional  uncertainty  into  those  calcula¬ 
tions  using  the  less  detailed  structural  models,  since  it  represents 
a  loss  of  information  calculated  with  FSI.  However,  it  is  crucial  if 
a  given  fire  scenario  is  to  be  associated  with  structural  collapse. 
These  issues  will  be  discussed  in  further  detail  below. 
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Fig.  4.  Portion  of  WTC  floor  gridded  at  resolution  used  for  FDS  simulations.  A  section  of  the  floor  truss  system  and  a  perimeter  column  cross-section  as  modeled  in  ANSYS 
are  also  shown.  The  truss  rods  are  2.5  cm  in  diameter  and  the  column  is  35  cm  on  a  side. 

Courtesy  K.  Prasad  (Baum,  2005). 


3.  Coupling  and  feedback 

The  sections  above  outline  how  the  overall  analysis  was  bro¬ 
ken  down  and  indicate  the  methods  used  to  study  each  portion  of 
the  overall  simulation  as  an  isolated  entity.  However,  if  the  effects 
of  a  given  fire  on  a  structure  are  to  be  predicted,  it  is  necessary  to 
couple  these  codes  in  an  internally  consistent  manner.  Since  the  air¬ 
craft  impact  is  both  peculiar  to  the  WTC  attack  and  considers  less 
than  1  s  of  elapsed  time,  attention  is  focused  on  the  interactions 
between  the  fire  dynamics,  heat  transfer,  and  structural  analysis. 
These  processes  occur  simultaneously  in  reality,  but  are  consid¬ 
ered  sequentially  in  the  WTC  analyses.  This  is  a  reflection  of  the 
current  state  of  the  art.  I  can  find  no  published  example  of  a  fully 
coupled  simulation  of  a  fire  in  a  complex  structure  that  includes 
fire  dynamics,  heat  transfer,  and  structural  mechanics. 

At  many  points  in  the  WTC  analyses,  photographic  and  visual 
evidence  are  used  to  compensate  for  our  current  inability  to  per¬ 
form  a  fully  coupled  simulation.  The  most  obvious  of  these  is  the 
effect  of  window  breaking  on  fire  dynamics.  While  glass  windows 
are  not  structural  components,  the  “collapse  analysis”  of  such  a 
window  requires  a  coupling  of  fire  dynamics,  heat  transfer  and 
stress  analysis  analogous  to  that  needed  for  the  load  bearing  struc¬ 
ture.  Fortunately,  there  is  no  need  to  consider  this  problem  in  detail 
since  it  was  the  subject  of  a  superb  recent  review  by  Pagni  (2003). 
From  the  present  perspective,  the  most  important  conclusion  to 
emerge  from  this  review  is  the  following:  While  the  time  to  initial 
crack  formation  in  the  glass  can  be  reliably  estimated  as  a  func¬ 
tion  of  window  geometry,  materials,  and  thermal  load,  the  time  at 
which  the  glass  falls  out  of  its  frame  cannot  be  predicted  with  any 
confidence.  Since  the  fire  ventilation  will  not  change  until  the  win¬ 
dow  pane  is  removed  from  its  frame,  this  aspect  of  the  feedback 
from  changes  in  the  building  structure  to  the  fire  can  not  yet  be 
quantified. 

3.2.  The  fire-structure  interface 

The  coupling  between  the  fire  dynamics  and  the  structure  is 
dominated  by  the  radiation  heat  transfer.  The  radiation  field  must 


be  determined  from  solutions  of  the  radiative  transport  equa¬ 
tion,  which  relates  the  incident  flux  to  the  spatial  distribution  of 
temperature  and  combustion  products  (most  particularly  the  dis¬ 
tribution  of  soot  particulate)  as  well  as  the  enclosure  geometry. 
Such  calculations  are  typically  performed  as  part  of  a  CFD  based 
simulation  of  the  fire  dynamics.  However,  the  ability  to  couple  such 
codes  as  the  NIST  Fire  Dynamics  Simulator  (FDS)  McGrattan  (2004) 
directly  to  a  suitable  structural  analysis  code  does  not  yet  exist.  The 
enormous  differences  in  spatial  and  temporal  length  scales,  differ¬ 
ences  in  numerical  techniques,  and  the  complexity  of  the  computer 
codes  makes  the  development  of  an  efficient  coupled  analysis  of 
fire-structure  interactions  a  daunting  task.  This  disparity  in  length 
scales  is  illustrated  in  Fig.  4,  which  shows  a  small  portion  of  a  WTC 
floor  gridded  at  the  50  cm  resolution  used  for  the  FDS  simulations 
together  with  a  small  segment  of  the  floor  truss  system  and  an  insu¬ 
lated  perimeter  column  as  represented  in  ANSYS.  Since  the  truss 
rods  have  a  2.5-cm  diameter  and  the  perimeter  columns  are  35  cm 
on  a  side  it  is  clear  that  FDS  must  be  supplemented  by  a  method 
that  permits  the  heat  transfer  to  the  structure  to  be  calculated. 

The  FSI  offers  an  approach  to  the  calculation  of  the  coupled  heat 
transfer  problem.  It  takes  advantage  of  the  fact  that  the  simplest 
compartment  fire  models,  the  “zone  models”,  assume  that  the  com¬ 
partment  is  divided  into  a  hot,  soot  laden  upper  layer  and  a  cool, 
relatively  clear  lower  layer.  For  the  purposes  of  the  FSI,  the  tem¬ 
perature  gradients  in  the  horizontal  directions  are  much  smaller 
than  those  in  the  vertical  direction.  The  layer  temperatures  and 
radiative  properties  are  taken  from  suitably  chosen  temporal  and 
spatial  averages  of  output  generated  by  FDS.  The  time  averages  are 
chosen  to  be  compatible  with  the  time  scales  associated  with  ther¬ 
mal  diffusion  through  the  smallest  structural  members  of  interest. 
The  spatial  averages  replace  the  detailed  vertical  temperature  and 
absorption  coefficient  profiles  with  an  effective  “zone  model”  pro¬ 
file.  This  last  approximation  is  not  crucial  to  the  development  of  the 
FSI,  but  is  used  to  reduce  the  information  transfer  required  between 
FDS  and  ANSYS. 

The  step  that  is  crucial  is  recognition  of  the  implications  of  the 
layered  thermal  and  optical  properties  of  the  gas  phase.  Let  length 
scales  in  the  horizontal  directions  x  be  scaled  with  a  length  L  (over 
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Fig.  5.  Predicted  upper  layer  temperature  of  a  floor  of  the  north  tower  30  min  after  impact  using  the  model  shown  in  Fig.  1  .  The  peak  temperatures  are  in  the  vicinity  of 
1100°C. 

Courtesy  McGrattan  (2004) 


60  m  in  the  WTC  towers),  while  the  grey  gas  absorption  coefficient 
k  and  length  scales  in  the  vertical  direction  z  are  scaled  with  H, 
the  height  of  an  individual  floor  (less  than  4  m  in  the  WTC  towers). 
Finally,  let  the  integrated  intensity  I(x ,  z,  fix,  fiz)  be  normalized 
with  respect  to  oT '4,  where  a  is  the  Stefan-Boltzmann  constant, 
fi  =  (fix,  fiz)  denotes  the  local  direction  of  the  radiation  field,  and 
Tr  a  suitable  reference  temperature.  Then,  denoting  dimensionless 
quantities  with  a  tilde,  the  radiative  transport  equation  takes  the 
form: 


H  -  „  ~  _  dl  „ 

—(fix  •  Vx7)  +  —  k 


Since  the  ratio  H/L  <<c  1 ,  the  terms  in  round  brackets  on  the  left  hand 
side  of  Eq.  ( 1 )  can  be  ignored.  The  remaining  terms  are  those  associ¬ 
ated  with  the  problem  of  radiative  transport  between  plane  parallel 
layers.  For  this  simplified  geometry,  the  radiative  transport  equa¬ 
tion  can  be  solved  exactly  and  explicit  formulae  for  the  heat  flux 
obtained  as  functions  of  the  temperatures,  hot  layer  depth,  soot 
concentration,  as  well  as  the  location  and  orientation  of  the  struc¬ 
tural  element.  The  basic  analysis  is  quite  well  known  in  the  heat 
transfer  community  (Siegel  and  Howell,  1992),  and  the  extensions 
needed  for  use  in  FSI  are  readily  obtained.  Details  can  be  found  in 
Prasad  and  Baum  (2005). 

The  FSI  converts  these  formulae  into  input  for  ANSYS  or  any 
other  finite  element  computer  code  that  can  be  used  for  both  time 
dependent  three  dimensional  heat  transfer  within  the  solid  phase 
and  for  the  stress  analysis  of  the  heated  structure.  The  solid  phase 
heat  transfer  is  dominated  by  thermal  conduction  with  thermal 
radiation  boundary  conditions.  The  input  radiation  is  calculated  by 
FSI  based  on  averages  extracted  from  FDS,  while  the  outgoing  radi¬ 
ation  is  a  function  of  the  surface  temperatures  predicted  by  ANSYS. 
Avery  detailed  geometry  model  of  all  solids  can  be  accommodated, 
with  solutions  that  account  for  insulation  damage  as  well  as  statis¬ 
tical  variations  in  the  application  of  spray-on  fireproofing.  Finally, 
the  FSI  interpolates  the  temperatures  obtained  from  this  calcula¬ 
tion  into  the  geometry  model  employed  for  the  stress  analysis.  The 
complexity  of  the  structural  mechanics  simulations  makes  this  last 
step  necessary. 


3.2.  Sample  simulations 

As  noted  above,  the  simulations  are  performed  sequentially, 
with  the  FDS  simulation  coming  first.  Fig.  5  shows  a  snapshot  of  the 
upper  layer  temperatures  in  the  north  tower  30  min  after  impact 
(McGrattan,  2004).  This  was  one  of  a  series  of  development  sim¬ 
ulations,  and  is  not  necessarily  the  best  prediction  of  the  state  of 
the  fires  at  that  time.  However,  it  gives  a  general  idea  of  the  kind 
of  predictions  that  can  be  made  using  this  methodology.  Note  that 
the  fire  dynamics  on  all  floors  must  be  computed  simultaneously, 
since  the  floors  are  connected  by  both  ventilation  shafts  and  impact 
damage.  The  results  show  fires  that  slowly  migrate  around  the  floor, 
with  peak  temperatures  in  the  upper  layer  in  the  vicinity  of  1 1 00  °C. 
Fig.  5  also  shows  the  relatively  smooth  temperature  variation  over 
most  of  the  upper  layer  with  rapid  transitions  confined  to  the  edges 
of  the  high  temperature  zones  where  burning  is  taking  place.  This 
helps  justify  the  plane  layer  radiative  transport  analysis  used  by  the 
FSI. 

The  computational  mesh,  even  at  50  cm  resolution,  requires 
nearly  150,000  cells  per  floor.  Over  500,000  time  steps  are  needed 
to  advance  the  calculation  for  6000  s  of  simulated  time.  Without 
the  use  of  the  parallel  processing  technique  incorporated  into  FDS 
as  part  of  the  WTC  investigation,  these  calculations  would  be  nearly 
impossible  to  perform.  Indeed,  McGrattan  (2004)  estimate  that  the 
computation  illustrated  in  Fig.  5  would  take  up  to  2  months  on  a 
single  processor,  if  it  were  capable  of  addressing  the  6-12  gigabytes 
of  memory  required  to  perform  the  computation.  This  memory 
requirement  would  rule  out  the  use  of  all  32  bit  “commodity”  pro¬ 
cessors.  The  use  of  parallel  processing  has  reduced  the  time  needed 
to  complete  a  simulation  to  as  little  as  two  days  using  eight  proces¬ 
sors  per  floor.  Moreover,  by  breaking  up  the  calculation  so  that  the 
computational  demands  on  each  processor  are  not  too  large,  inex¬ 
pensive  commodity  processors  in  “Beowulf  Clusters”  connected 
with  gigabit  networks  can  be  employed. 

Since  the  time  step  used  by  FDS  is  of  order  0.01  s,  and  the 
smallest  time  scales  of  interest  in  the  thermal  analysis  are  several 
seconds,  only  averaged  values  of  T4  and  the  absorption  coefficient 
at  several  second  intervals  need  be  passed  from  FSI  to  ANSYS.  It 
also  opens  the  possibility  of  coupling  the  fire  and  thermal  analy¬ 
ses  together.  So  long  as  hundreds  of  FDS  updates  can  be  performed 
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Fig.  6.  Predicted  temperatures  in  the  steel  structure  of  the  96th  floor  of  the  north  tower  100  min  after  impact.  The  peak  temperatures  reach  675  °C.  The  calculation  for  each 
floor  includes  the  insulation  and  concrete  slab  (not  shown). 

Courtesy  K.  Prasad  (Baum,  2005). 


between  solid  phase  thermal  updates,  the  cost  of  restarting  FDS 
with  a  corrected  set  of  thermal  boundary  conditions  derived  from 
the  solid  phase  analysis  can  be  minimized.  This  issue  will  be 
addressed  again  below. 

Predicted  temperatures  in  the  steel  structural  components  of 
the  96th  floor  of  the  north  tower  1 00  min  after  impact  are  shown  in 
Fig.  6.  This  is  one  of  a  series  of  development  simulations,  and  is  not 
necessarily  the  best  prediction  of  either  the  thermal  or  mechanical 
state  of  the  structure  at  that  time.  However,  it  does  illustrate  the 
extent  to  which  the  FSI  permits  highly  detailed  calculations  of  the 
thermal  state  of  the  structure  to  be  coupled  with  FDS  simulations. 
The  peak  temperatures  in  this  simulation  reach  675°  C  near  the 
south  face  of  the  floor,  opposite  the  damaged  portion  of  the  floor. 
This  is  roughly  the  area  reached  by  the  fires  several  minutes  earlier 
as  they  moved  about  the  floor.  The  simulation  includes  the  concrete 
slab  and  the  insulation  material,  which  are  omitted  for  clarity. 

The  above  simulation  spans  a  range  of  length  scales  ranging  from 
over  60  m  (the  tower  width)  to  1  cm  (the  radius  of  a  truss  rod).  The 
entire  floor  system  of  each  individual  floor,  including  insulation 
and  the  concrete  floor  slab,  is  coupled  together  in  a  single  calcula¬ 
tion.  The  ANSYS  model  for  a  single  combined  floor  system  contains 
approximately  450,000  nodes  and  650,000  elements  (Prasad  and 
Baum,  2005).  Each  column  is  simulated  separately,  since  they  are 
weakly  coupled  thermally.  The  column  models  have  20  nodes  per 
floor  in  the  vertical  direction  and  a  few  hundred  to  describe  the 
heat  transfer  in  each  horizontal  cross  section.  Thus,  the  number  of 
nodes  and  elements  in  the  columns  is  comparable  to  that  in  the 
floor  system.  It  takes  approximately  one  day  per  floor  to  perform 
a  thermal  analysis.  Clearly,  calculations  at  this  level  of  detail  are 
only  possible  because  of  the  relative  simplicity  of  conduction  heat 
transfer  compared  with  turbulent  combustion.  The  fact  that  each 
floor  can  be  simulated  independently  is  itself  a  form  of  parallel  pro¬ 
cessing.  It  requires  many  copies  of  the  simulation  code  (ANSYS  in 
this  case),  as  opposed  to  a  single  copy  of  FDS  for  the  fire  simula¬ 
tions.  However,  financial  issues  aside,  it  leads  to  highly  efficient 
computations. 

The  efficiency  of  the  computations  is  of  more  than  passing  inter¬ 
est.  Given  the  propagation  of  uncertainty  through  the  chain  of 
simulations,  it  is  not  possible  to  obtain  credible  results  without 


many  repetitions  of  each  calculation.  Only  in  this  way  can  the 
sensitivity  of  the  calculations  to  plausible  variations  in  input  be 
established.  Moreover,  the  only  way  to  make  contact  with  the  pho¬ 
tographs  and  videos  that  provide  a  check  on  the  credibility  of  the 
simulations  is  to  repeat  them  enough  times  to  establish  the  plau¬ 
sibility  of  the  results.  Tens  of  calculations  of  the  type  illustrated  in 
the  figures  shown  here  have  been  performed  as  part  of  the  NIST 
WTC  investigation. 

One  novel  feature  of  the  thermal  analysis  performed  by  the  FSI  is 
worth  special  mention.  The  thickness  of  the  spray-on  insulation  on 
the  trusses  in  the  WTC  towers  has  been  the  focus  of  considerable 
attention  (Quintiere  et  al.,  2002).  However,  given  the  method  of 
application,  the  thickness  can  only  be  defined  statistically,  as  there 
is  a  wide  variation  in  the  actual  thickness  at  each  point  along  the 
protected  components.  Fig.  7  illustrates  the  methodology  devel¬ 
oped  to  deal  with  this  fact.  An  artificially  thick  uniform  insulation 
model  is  built  up  around  each  protected  rod.  Then  a  mean  thick¬ 
ness  and  variance  are  chosen,  and  individual  elements  are  chosen 
statistically  to  have  either  the  proper  thermal  properties  or  essen¬ 
tially  no  thermal  inertia.  The  result  is  a  finite  element  model  whose 
geometric  simplicity  is  preserved  but  which  accounts  for  the  loss  of 
protection  associated  with  random  thickness  variations.  The  result 
is  a  net  loss  for  a  given  average  thickness  because  those  portions 
of  the  rod  with  less  protection  experience  higher  heat  transfer.  The 
excess  heat  can  be  transferred  a  small  distance  along  the  rod  by 
axial  diffusion,  canceling  whatever  reduction  might  have  resulted 
from  the  thicker  adjacent  insulation.  This  also  gives  rise  to  local 
“hot  spots”,  which  can  cause  a  weakening  of  the  structural  element 
in  the  immediate  vicinity  of  the  spot. 

4.  Research  issues 

The  remaining  issues  to  be  discussed  here  concern  the  coupling 
between  the  thermal  analysis  of  the  structure  and  the  calculation 
of  stress  and  displacement  fields.  While  they  may  not  be  “research 
issues”  in  the  strictest  sense  of  the  word,  the  extent  to  which  these 
phenomena  need  to  be  coupled  is  not  adequately  appreciated  in 
the  view  of  this  author.  The  first  point  worth  emphasizing  is  that 
the  spatial  and  temporal  resolution  of  the  stress  analysis  cannot  be 
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Fig.  7.  Finite  element  model  of  spray-on  insulation  showing  random  variations  in  thickness  (top).  The  outer  layer  of  cells  have  no  thermal  inertia,  the  middle  layer  represents 
the  actual  insulation  whose  thickness  is  specified  statistically,  and  the  inner  layer  simulates  the  steel  rod.  The  lower  figure  shows  the  computed  temperatures. 

Courtesy  K.  Prasad  (Baum,  2005). 


Fig.  8.  Kernel  functions  normalized  by  their  value  at  ri  =  0.  The  quantities  actually 
plotted  are  G(r])  =  G(rj,  r)/G(0,  t)  and  K(ri)  =  K(rj,  r)/I<(0,  r).The  Gaussian  function 
G  decays  more  rapidly  than  K  which  decays  ~?7-3. 


chosen  arbitrarily  once  the  requirements  for  thermal  analysis  have 
been  decided. 

Consider  the  idealized  problem  of  a  three-dimensional  elastic 
half-space  that  is  loaded  thermally  by  a  moving  point  source  of  heat 
q0Qj{r)  concentrated  at  a  point  on  the  free  surface  r  =  R( r).  Then, 
the  classic  theory  of  linear  thermo-elasticity  can  be  used  to  show 
that  without  approximation  the  surface  temperature  rise  and  the 
normal  component  of  the  surface  displacement  (the  “bulge”)  take 
the  following  form  (see  Fig.  8): 


0(r,  r) 


dr0Qj(T0)G(T 


t o , 


r  -  R(t0) 

\Ar-ro) 


w(r,  t) 


dr0Qj{r0)K{r 


r0, 


r-R{  r0) 
^{r-r0) 


Here,  G  is  the  classical  Gaussian  Kernel  function  and  I<  is  an  anal¬ 
ogous  function  that  emerges  from  the  thermo-elastic  analysis.  The 
variable  r  is  a  time  normalized  with  the  time  scale  t0  imposed  by 
the  thermal  heating,  and  the  length  scale  /  =  y7 t0/a  where  a  is 
the  thermal  diffusivity.  The  dimensionless  temperature  rise  0  is 


made  non-dimensional  with  respect  to  q0llk  where  k  is  the  ther¬ 
mal  conductivity.  The  surface  displacement  w  is  normalized  by  l 
and  a  thermo-elastic  “coupling  constant”.  The  function  I<  can  be 
written  explicitly  as: 

*--(¥£)  3^  **(£)  («•(?)-'■(?))  <«> 

The  quantities  X  and  /x  are  the  Lame  constants  of  classical  elastic¬ 
ity,  while  J0  and  l\  are  Modified  Bessel  functions.  These  results  are 
probably  well  known,  although  I  cannot  find  them  in  the  classical 
literature.  A  derivation  of  these  results  is  presented  in  the  appendix. 

The  importance  of  these  formulae  lies  in  the  fact  that  both  the 
temperature  rise  and  displacement  are  shown  to  be  inherently 
time-dependent  diffusive  processes,  with  the  time  scale  and  dif¬ 
fusion  set  by  the  temperature  field.  Moreover,  the  displacements 
depend  explicitly  on  the  entire  space-time  history  of  the  thermal 
heating  of  the  exposed  surface.  If  the  temperature  field  is  indepen¬ 
dent  of  time  or  if  it  is  postulated  as  opposed  to  derived,  then  this 
issue  does  not  arise.  However,  under  either  of  these  circumstances, 
it  is  no  longer  possible  to  couple  the  fire  dynamics  to  the  structural 
analysis.  While  the  problems  of  interest  in  a  structural  collapse  sce¬ 
nario  are  inherently  non-linear,  that  does  not  invalidate  the  above 
conclusions.  The  thermally  induced  diffusion  is  then  superimposed 
on  a  variety  of  other  phenomena. 

The  significance  of  geometric  non-linearity  caused  by  large  ther¬ 
mally  induced  displacements  has  recently  been  emphasized  by 
Lane  (2003).  This  has  important  consequences  for  coupled  thermal- 
structural  analyses.  Even  if  large  deflections  in  components  of  the 
structure  do  not  change  the  fire  dynamics  appreciably,  the  defor¬ 
mation  can  change  the  thermal  environment  experienced  by  the 
component.  Fig.  9  shows  a  comparison  between  two  time  depen¬ 
dent  ANSYS  simulations  of  a  horizontal  rod  fixed  at  the  ends  and 
placed  in  a  gas  with  a  temperature  distribution  that  increases  with 
increasing  height.  The  simulation  on  the  left  assumes  the  heat 
transfer  to  the  rod  is  evaluated  at  its  original  position  while  the 
right  hand  simulation  uses  the  current  position  of  the  rod.  The 


Fig.  9.  Temperatures  (top)  and  displacements  of  an  unprotected  steel  rod  in  a  gas  with  a  vertically  stratified  temperature  profile.  The  simulation  on  the  left  assumes  the  heat 
transfer  to  the  rod  is  evaluated  at  its  original  position  while  the  right  hand  simulation  uses  the  current  position  of  the  rod. 

Courtesy  K.  Prasad  (Baum,  2005). 
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upper  part  of  the  figure  shows  the  temperatures  while  the  deflec¬ 
tions  are  shown  below.  Note  that  the  temperature  in  the  left  hand 
simulation  is  uniform,  because  the  heat  flux  to  the  rod  is  always 
evaluated  at  the  original  location,  even  though  the  temperature 
in  the  rod  increases  with  time.  The  change  in  heat  flux  to  the  rod 
with  changing  location  is  included  in  the  right  hand  simulation. 
As  a  result,  the  rod  temperature  varies  significantly  with  position. 
The  deflections  in  this  case  are  smaller,  because  as  the  rod  deforms 
it  moves  into  a  cooler  environment  with  lower  heat  transfer  as  a 
result.  While  this  may  not  seem  too  important  for  an  isolated  truss 
rod,  the  consequences  if  repeated  on  a  large  scale  can  be  quite  sig¬ 
nificant  (Usmani  etal.,  2003;  Lane,  2003).  Both  effects  illustrated  in 
this  section  could  be  addressed  if  thermal  and  stress  analyses  of  a 
given  structure  were  performed  in  parallel  rather  than  sequentially. 
The  difficulty  arises  in  choosing  a  consistent  level  of  modeling  detail 
which  preserves  the  ability  to  perform  the  stress  analysis  without 
compromising  the  accuracy  of  the  temperature  field. 


that  this  continuous  change  must  affect  the  stresses  is  ignored  in 
this  approach.  This  technique  is  justified  by  noting  that  the  elastic 
wave  propagation  speed  is  so  fast  compared  with  the  time  scales 
of  interest  in  thermo-elastic  phenomena  induced  by  fires  that  a 
quasi-steady  analysis  is  justified. 

The  analyses  that  follow  are  intended  to  show  that  computa¬ 
tional  techniques  that  “freeze”  the  temperature  at  a  given  time 
and  compute  an  equilibrium  stress  distribution  are  not  consis¬ 
tent  with  the  dynamical  equations  of  thermo-elasticity,  even  if  the 
elastic  wave  propagation  speed  is  taken  to  be  infinite.  The  next  sec¬ 
tion  demonstrates  how  the  general  solutions  to  the  equations  of 
thermo-elasticity  couple  the  time  scale  for  the  evolution  of  the  dis¬ 
placements  to  that  of  the  temperature  field.  In  particular,  it  is  shown 
that  the  solutions  for  the  displacements  cannot  obey  an  equilibrium 
equation  unless  the  temperature  field  is  independent  of  time.  Fol¬ 
lowing  this,  formal  solutions  for  a  half-space  loaded  thermally  are 
derived.  Again,  it  is  clear  that  part  of  the  solution  for  the  stresses 
and  displacements  are  inherently  time  dependent. 


5.  Concluding  remarks 

A  discussion  of  the  current  state  of  our  ability  to  predict  the 
effects  of  fire  on  building  structures  has  been  presented.  The 
emphasis  has  been  on  the  degree  to  which  fire  dynamics,  heat 
transfer,  and  structural  analysis  can  be  coupled.  The  NIST  inves¬ 
tigation  into  the  collapse  of  the  WTC  towers  is  used  as  a  backdrop, 
because  it  provides  a  concrete  illustration  of  the  strengths  and  limi¬ 
tations  of  our  current  methodology.  While  the  simulation  capability 
in  each  of  these  fields  can  (and  no  doubt  will)  be  advanced  sepa¬ 
rately,  it  is  only  when  they  are  coupled  together  in  some  way  that 
fire  effects  on  structures  can  be  quantitatively  assessed.  This  is  an 
area  that  would  benefit  greatly  from  increased  attention  by  the 
Fire  Research  community.  I  hope  this  community  will  rise  to  the 
challenge. 


Appendix  A.  Thermal  stress  analysis 

AA.  Background 

The  following  analysis  is  motivated  by  the  observation  that 
the  calculation  of  stresses  induced  by  time  dependent  tempera¬ 
ture  fields  using  certain  commercial  software  packages  seems  to 
be  much  more  difficult  than  the  corresponding  calculation  for  an 
isothermal  environment.  While  the  applications  of  interest  clearly 
involve  highly  non-linear  calculations,  the  starting  point  for  most 
fire  scenarios  is  almost  always  an  undamaged  building  at  room 
temperature.  Since  virtually  all  buildings  are  designed  to  keep  the 
stresses  well  below  the  elastic  limit  and  the  deflections  of  the  load 
bearing  structure  reasonably  small,  the  starting  point  for  simula¬ 
tions  of  fire  induced  damage  must  lie  within  the  domain  of  linear 
elasticity.  Moreover,  the  difficulties  that  arise  using  the  packages  of 
interest  are  evident  before  the  temperature  rise  is  large  enough  to 
affect  the  elastic  or  thermal  properties  of  most  structural  materi¬ 
als.  Under  these  circumstances  the  thermally  induced  stresses  are 
also  linear,  and  the  temperature  fields  can  be  described  by  the  heat 
conduction  equation  for  the  material(s)  of  interest. 

The  facts  described  above  justify  an  analysis  of  the  coupling 
between  the  temperature  and  thermally  induced  stresses  based  on 
the  linear  thermo-elastic  equations.  The  temporal  dependence  of 
the  stresses  is  the  focus  of  this  analysis.  A  popular  technique  for 
solving  the  thermo-elastic  equations  is  to  first  compute  (or  assume) 
the  time  dependent  temperature  distribution  in  the  load  bearing 
structure.  Then,  given  this  information,  the  temperature  distribu¬ 
tions  are  “frozen”  at  a  succession  of  discretely  chosen  times  and  an 
equilibrium  solution  is  sought  for  the  state  of  stress  at  each  chosen 
time.  The  fact  that  the  temperature  is  changing  continuously  and 


A.2.  The  thermo-elastic  equations 


The  starting  point  for  the  analysis  are  the  linear  thermo-elastic 
equations  which  relate  the  displacements  ut  and  stresses  r,j  to  each 
other  and  to  the  temperature  T  in  the  elastic  medium.  They  can  be 
written  in  the  form  Sokolnikoff  (1956): 


=  Fi  + 


d rij 
dxj 


(A.1) 


a(3X  +  2/n)(T  -T0)8ij 


(A.2) 


Here,  p  is  the  density  of  the  material,  F,-  is  the  body  force  per  unit 
volume,  pt  and  X  the  Lame  constants,  a  the  coefficient  of  thermal 
expansion,  and  T  the  temperature.  The  temperature  T0  is  the  ref¬ 
erence  temperature  of  the  material  in  its  unstressed  state  before 
the  fire.  That  temperature  is  taken  to  be  uniform  here,  which  is 
a  reasonable  simplification  given  the  temperature  rise  associated 
with  a  building  fire.  These  equations  are  supplemented  by  suitable 
boundary  conditions  that  express  the  connections  to  other  portions 
of  the  structure  and  external  loads.  The  temperature  field  evolves 
according  to  the  heat  conduction  equation. 


d(T-T0) 

3t 


d2(J -T0) 

K 


(A.3) 


The  coupling  between  the  stresses  and  the  temperature  field 
comes  from  the  effect  of  the  temperature  gradients  on  the  volu¬ 
metric  expansion  0.  The  body  force  plays  no  role  in  this  and  will 
henceforth  be  ignored.  Taking  the  divergence  of  the  thermo-elastic 
evolution  Eq.  (A.l )  yields: 


(A.4) 


(X  +  2 n)V2(j>  -  a( 3X  +  2/i.)V2(T  -  T0) 


(A.5) 


The  next  step  is  to  simplify  these  equations  by  eliminating 
the  fast  wave  motion  associated  with  the  irrotational  waves.  It 
is  easy  to  see  from  Eq.  (A.5)  that  these  waves  propagate  with  a 
speed  c2  =  {X  +  2pi)l p.  For  typical  steels,  c~0(103)m/s.  Moreover, 
the  thermal  diffusivity  k/(pCp)  ~  0(10-5)  m2/s.  We  now  introduce 
dimensionless  thermal  and  mechanics  variables  respectively  as  fol¬ 
lows: 


T  -T0 


®{yi,  r). 


=  (y„  t 


(A.6) 


H.R.  Baum  /  Mechanics  Research  Communications  38  (201 1)  1-11 


9 


*-**>,,+  («) 

(1+2/0  (A'8) 

The  non-dimensional  variables  are  chosen  so  that  the  full  time 
dependent  heat  conduction  equation  is  retained,  and  the  temper¬ 
ature  rise  is  related  to  a  heat  flux  to  a  bounding  surface  q0.  The 
coupling  parameter  /3  scales  the  thermally  induced  deformations 
to  the  temperature  rise,  and  e  is  the  ratio  of  the  speed  of  a  ther¬ 
mal  “front”  to  the  irrotational  wave  speed  in  the  elastic  medium. 
The  time  scale  t0  is  arbitrary,  and  its  choice  sets  both  the  diffusion 
controlled  length  scale  and  the  magnitude  of  e.  Given  the  approxi¬ 
mate  values  of  wave  speed  and  thermal  diffusivity  shown  above,  it 
is  clear  that  6  <1  for  any  time  scale  of  interest  in  a  fire  scenario. 

The  dimensionless  evolution  equations  for  the  normalized  dila¬ 
tion  0  and  temperature  rise  0  take  the  following  form: 


90  =  V2© 

9r 

(A.9) 

9  920  9  9 

e\?=Vcp  V0 

9r2 

(A.10) 

> 

II 

e 

(A.  11) 

It  is  now  easy  to  show  the  structure  of  the  solution  for  the  irrota¬ 
tional  portion  of  the  deformation.  Ignoring  terms  of  order  e2,  it  is 
clear  that  0  must  have  the  form: 

0  =  0  +  0*,  V20*=O  (A.12) 


must  miss  at  least  part  of  the  solution  to  the  equations  of  thermo¬ 
elasticity.  Eq.  (A.15)  shows  that  when  the  temperature  is  “frozen”, 
the  right  hand  side  of  the  equation  is  replaced  by  zero  over  almost 
all  time,  except  at  the  discrete  points  where  a  Dirac  delta  func¬ 
tion  is  in  effect  used  to  update  the  temperature  field.  The  accuracy 
of  this  kind  of  approximation  is  completely  unknown.  Eq.  (A.17) 
shows  that  the  non-equilibrium  portion  of  the  solution  evolves  on 
the  same  time  scale  as  the  solution  to  the  heat  conduction  equa¬ 
tion.  Moreover,  since  the  heat  conduction  solutions  are  smooth, 
the  solution  to  the  non-equilibrium  displacement  field  must  also 
be  smooth.  This  would  seem  to  rule  out  approximating  the  time 
dependence  of  the  temperature  field  by  a  discrete  set  of  steps  for 
the  purpose  of  calculating  the  displacement  fields. 

A3.  Heated  elastic  half-space 

In  order  to  make  some  of  these  ideas  more  precise,  consider  the 
idealized  problem  of  an  elastic  half-space  heated  at  the  surface  by 
a  prescribed  heat  flux  that  in  general  depends  on  space  and  time. 
The  variables  are  made  dimensionless  as  before,  with  the  coordi¬ 
nate  normal  to  the  surface  pointing  into  the  solid  denoted  by  z,  and 
the  coordinates  parallel  to  the  surface  denoted  by  r  =  (x,y).  The 
dimensional  heat  flux  to  the  surface  at  z  =  0  is  given  by  qz  =  q0 Q(x ,  y, 
r).  Furthermore,  body  forces  are  ignored  and  there  are  no  mechan¬ 
ical  forces  acting  on  the  surface  z  =  0.  Thus,  the  only  reason  any 
stresses  and  deformations  are  set  up  in  the  solid  is  because  of  the 
thermal  loading. 

Under  these  circumstances,  the  boundary  conditions  at  the  sur¬ 
face  become: 


This  clearly  shows  that  the  dilation  has  two  parts:  a  harmonic  con¬ 
tribution  0*  that  can  be  regarded  as  a  local  equilibrium  solution 
corresponding  to  the  instantaneous  boundary  conditions,  and  a 
part  that  is  directly  proportional  to  the  local  temperature  rise.  This 
part  of  the  solution  never  is  in  equilibrium,  unless  the  temperature 
is  in  steady  state.  While  a  steady  temperature  field  may  be  a  useful 
approximation  to  the  state  of  a  furnace,  for  example,  it  is  a  very 
poor  representation  of  the  temperature  distribution  resulting  from 
a  building  fire. 

This  result  can  be  seen  even  more  clearly  by  noting  that  since 
the  dimensionless  displacement  is  a  vector  field,  it  can  always  be 
decomposed  into  an  irrotational  and  a  solenoidal  part. 

v  =  W  +  V  x  A  (A.13) 

Then,  since  V  •  v  =  0  =  V2^,  the  scalar  potential  function  0  satis¬ 
fies  the  equation: 

V20  =  0  +  0*  (A.14) 

The  harmonic  function  0 *  can  be  eliminated  to  obtain  an  explicit 
relation  between  the  irrotational  component  of  the  displacement 
and  the  temperature  field  in  the  following  form: 

9  9  90 

SJ2W20=—~  (A.15) 

dr 

Here,  the  fact  that  0  satisfies  the  heat  conduction  equation  has  been 
used.  Alternatively,  the  solution  for  0  can  be  decomposed  into  an 
evolution  equation  and  an  equilibrium  equation  as  follows: 


0  =  0*  +  ^ 


(A.16) 


V2V20*  =  0, 


(A.17) 


Both  Eqs.  (A.15)  and  (A.17)  show  that  only  part  of  the  solution 
for  the  thermally  induced  displacement  can  correspond  to  a  local 
equilibrium  state  if  the  temperature  varies  with  time.  Any  solution 
procedure  that  works  by  “freezing”  the  temperature  in  time  and 
using  an  equilibrium  equation  to  satisfy  the  boundary  conditions 


dry 

rzi  =  0,  i  =  x,y,  z;  —  =  -Q(x,y,  r)  (A.18) 

The  stress  free  boundary  conditions  can  be  conveniently  rewritten 
in  terms  of  displacements  by  introducing  the  parallel  displacement 
vector  V  =  ( u ,  v),  the  perpendicular  displacement  w,  and  the  paral¬ 
lel  gradient  operator  Wh  defined  by: 


(  d_  _9_ 
[dx ’  9y 


(A.19) 


The  requirement  that  the  two  shear  components  of  the  stress  van¬ 
ish  at  the  surface  can  then  be  written  in  the  form: 

Vh  w  +  —  V  =  0  @z  =  0  ( A.20) 

dz 

Similarly,  the  normal  stress  will  vanish  at  the  surface  provided  that: 

r\ 

X0  +  2ii^  =  {X  +  2/x)0  @z  =  0  (A.21 ) 

In  the  present  notation,  0  can  be  written  in  the  form: 

r\ 

0  =  Vh-\/+^  =  0  +  0*  (A.22) 

dz 

The  heat  flux  is  assumed  to  be  applied  to  the  surface  over  a  finite 
area.  Thus,  the  temperature  rise  together  with  all  displacements 
must  vanish  as  z^  oo  and  r  ->  oo. 

The  equations  that  must  be  solved  can  now  be  rewritten  as 
follows.  The  heat  conduction  equation  takes  the  form: 


90 

9r 


0 


(A.23) 


Using  the  general  results  obtained  in  the  previous  section,  the  equi¬ 
librium  equations  for  the  parallel  components  of  the  displacement 
become: 

(v^  +  £)'?+(^)v',0*  =  Vh0  (A-24) 
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The  auxiliary  function  0*  defined  in  Eq.  (A.12)  satisfies  the  Laplace 
equation,  written  in  the  present  notation  as: 


( 


V^  + 


0*  =  0 


(A.25) 


The  system  of  equations  that  must  be  solved  thus  consists  of  Eq. 
(A.23),  which  determines  0,  Eq.  (A.24),  which  determines  V ,  Eq. 
(A.25),  which  determines  0*,  and  Eq.  (A.22),  which  determines  w. 

The  solutions  can  be  obtained  using  transform  methods.  Let 
the  Fourier-Laplace  transform  of  an  arbitrary  function  /  (r,  z,  r) 
be  defined  as  follows: 

/oo  poo 

d2r  /  dr  exp(-ik  ■  r  -pr)/(r,  z,  r)  (A.26) 
oo  J  0 


While  this  is  not  a  realistic  representation  of  the  spatial  distribu¬ 
tion  of  the  heat  flux  induced  by  an  individual  fire  to  a  large  floor 
area,  it  does  pick  up  two  key  features  of  such  a  fire.  First,  the  overall 
strength  of  the  fire  (measured  by  its  overall  heat  release  rate)  will 
change  with  time.  Second,  the  fire  will  migrate  from  place  to  place 
as  its  fuel  is  consumed  and  the  availability  of  oxygen  changes  with 
time.  Under  these  circumstances  the  surface  temperature  distribu¬ 
tion  simplifies  to: 


@(r,  r,0)  =  /  dr0 - Qr(ro) 

'o  4(7r(r-r0))3/2 


exp 


m  = 


r-R(  To) 


\/0-To) 


(A.37) 

(A.38) 


The  Fourier-Laplace  transform  of  the  solutions  satisfying  boundary 
conditions  at  infinity  can  then  readily  be  found  to  be: 


0  =  Q{k,p) 


exp  (-\/P  +  fe2Z) 


y/p  +  k2 


0*  =  A(k,  p)  exp (-kz) 


V: 
0  = 

W= 

CO  = 


ik&(k,  p ) 

:  0  +  B(/c,p)exp(-/<z)  + 


(A.27) 


(A.28) 

(A.29) 


A(/c,p)—  exp(-kz)  (A.30) 


-  ^  exp  \Jp  +  /<2z^  -kB  exp  (-kz)-coA  exp(-kz)  (A.31) 
3  jjj  -|-  A,  -\~  ( /x  +  A.)  kz 


2  k/x 


(A.32) 


The  unknown  functions  A(/c,p)  and  B(k,p)  are  determined  by 
the  surface  boundary  conditions  given  in  Eqs.  (A.20)  and  (A.21). 
The  results  are: 


A  = 


li 


2  kQ, 


/x  +  A. 


1  - 


k 


B  = 


d 


(/X  +  A.)p  l  y^p  +  k2 


Cp  +  k2 

2/i  +  X 

k 


(A.33) 


(A.34) 


Since  the  primary  interest  in  this  solution  is  the  extent  to  which 
the  time  dependence  is  imbedded  in  the  result,  attention  is  focused 
on  the  vertical  displacement  at  the  surface.  This  part  of  the  solution 
can  be  readily  obtained  and  interpreted  in  the  light  of  the  general 
formulation  discussed  in  the  previous  section.  Physically,  it  repre¬ 
sents  the  “bulge”  in  the  surface  that  would  appear  in  the  vicinity 
of  the  heated  area.  The  recipe  for  the  bulge  can  be  readily  com¬ 
pared  with  that  for  the  temperature  rise  at  the  surface  induced 
by  the  heat  transfer.  Since  it  is  well  known  that  the  temperature 
distribution  must  be  treated  as  a  transient  phenomenon,  the  sim¬ 
ilarities  and  contrasts  between  these  two  results  will  lend  insight 
into  the  importance  of  a  coupled  transient  analysis  of  the  stresses 
and  displacements  induced  by  the  heat  transfer. 

First  consider  the  surface  temperature  distribution.  Using  the 
convolution  theorem,  the  solution  can  be  written  in  the  form: 


0(f,  T,  0) 


r 0)C(r-r0, 


(A.35) 


C(r,  r) 


1 

4(^r)3/2 


exp 


(A.36) 


In  order  to  simplify  the  analysis,  consider  the  special  case  where 
Q(r,  r)  is  concentrated  at  a  point  r  =  R( r)  with  a  strength  Qr(r). 


The  most  relevant  parts  of  this  result  in  the  present  context  are 
the  dependence  of  the  solution  on  the  previous  history  of  the  sur¬ 
face  heat  flux  distribution,  and  the  fact  that  the  Greens  function  G 
has  a  structure  determined  primarily  by  the  similarity  variable  r\j 
which  itself  is  inherently  time  dependent.  Note  that  the  singular¬ 
ity  in  the  integrand  at  r  =  ft(r0),  r  =  r0  is  an  artifact  of  injecting  a 
finite  heat  flux  into  a  point.  The  results  make  perfectly  good  sense 
for  points  away  from  the  current  location  of  the  heat  flux  source. 

The  solution  for  the  displacement  normal  to  the  surface  can  be 
found  in  an  analogous  manner.  However,  since  the  inversion  pro¬ 
cess  is  somewhat  more  complex  here  a  few  details  are  provided. 
The  Fourier-Laplace  transform  of  the  surface  displacement  takes 
the  form: 


w(k,p,0) 


2/i  +  l\  Q  |  k 
/x  +  A.  y  p  1  y/p  +  k2 


(A.39) 


The  solution  for  w  (r ,  r,  0)  takes  the  same  form  as  Eq.  (A.35)  except 
that  the  Greens  function  G  (r,  t)  is  replaced  by  a  new  kernel  func¬ 
tion  I<(r ,  r)  defined  as: 


I< 


/  2/x  +  A.  \ 
V  /x  T  A,  / 


d2k 

(2 Tj2 


dp  exp(pr) 
2  m  p 


(A.40) 


Carrying  out  the  inversion  of  the  Laplace  transform  first,  and  noting 
that  the  resulting  expression  depends  only  on  k  =  \k\,  the  Fourier 
inversion  integral  can  be  reduced  to: 

K  =  -(jrrr)^  l  ^erfc(?)j0^),  v  =  j=  (a.4d 

Here,  Jo  denotes  the  Bessel  function  of  the  first  kind  of  order  zero. 
The  final  integral  can  be  evaluated  with  the  aid  of  Mathematica  to 
yield: 


I<  = 


2/x  +  AA  1 

— - t —  — —  exp 

/x  +  X  J  8ttt 


(A.42) 


The  quantities  J0  and  Ji  are  the  modified  Bessel  functions  of  the  first 
kind  of  order  zero  and  one  respectively.  The  minus  sign  in  front  of 
the  (positive)  expression  for  I<  arises  from  the  definition  of  positive 
w  pointing  into  the  material.  Thus,  the  thermal  expansion  induces 
a  bulge  out  of  the  plane  of  the  surface,  so  that  w  must  be  negative. 
Finally,  if  the  heat  flux  to  the  surface  is  concentrated  at  a  point  as 
described  above,  the  solution  for  the  normal  surface  displacement 
becomes: 


w=  dToQr{ro)K(T-To,rjT ), 

Jo 


rjT 


f-R(T0) 
^/(r  —  r0) 


(A.43) 
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Clearly,  the  only  significant  differences  between  the  solutions 
for  &  and  w  at  the  surface  are  in  the  mathematical  structure  of  the 
kernel  functions  G  and  K.  They  are  plotted  in  Fig.  8  in  normalized 
form,  so  that  the  value  at  r]  =  0  for  each  function  is  one.  The  thermal 
function  G  has  an  exponential  decay,  consistent  with  the  inherently 
transient  diffusion  of  heat  from  the  point  source  at  the  surface.  The 
normal  displacement  kernel  I<  however,  decays  algebraically  with 
large  r\  with  This  is  a  consequence  of  the  fact  that  the  dis¬ 

placements  have  both  equilibrium  and  transient  components.  On 
the  other  hand,  the  spatial  dependence  of  both  functions  appears 
in  the  inherently  transient  independent  variable  rj  =  r/Vr,  which 
describes  a  diffusion  controlled  process.  Moreover,  the  full  solu¬ 
tions  for  both  the  temperature  and  the  displacement  are  dependent 
on  the  previous  history  of  each  quantity. 
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